Phase diagram, correlation gap, and critical properties of the Coulomb glass 
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We investigate the lattice Coulomb glass model in three dimensions via Monte Carlo simulations. 
No evidence for an equilibrium glass phase is found down to very low temperatures, although the 
correlation length increases rapidly near T = 0. A charge-ordered phase (COP) exists at low disorder. 
The transition to this phase is consistent with the Random Field Ising universality class, which shows 
that the interaction is effectively screened at moderate temperature. For large disorder, the single- 
particle density of states near the Coulomb gap satisfies the scaling relation g(e,T) — T s f(\e\/T) 
with S — 2.01 ± 0.05 in agreement with the prediction of Efros and Shklovskii. For decreasing 
disorder, a crossover to a larger effective exponent occurs due to the proximity of the COP. 

PACS numbers: 64.70ph,71.23.-k,75.10.Nr 
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In disordered insulators, the localized electrons cannot 
screen effectively the Coulomb interaction at low tem- 
perature. Therefore, many-electron correlations are im- 
portant in this regime. The long-range repulsion induces 
a soft "Coulomb gap" in the single-particle density of 
states (DOS). Efros and Shklovskii (ES) Q argued that 
the gap has a universal form g(e) oc |e — fi\ 5 near the 
chemical potential /i, with 5 > d— 1 in d dimensions, and 
that a saturated bound 5 = d — 1 modifies the variable- 
range hopping resistivity lni? ~ T~ x from Mott's law 
x = l/(d+ 1) to x = 1/2. Both the existence of the 
gap and the crossover to x ~ 1/2 at low temperature 
T have been confirmed experimentally and in numerical 
simulations Q , but the validity of 8 — 2 for d — 3 has yet 
to be firmly established. Pseudo ground-state numerical 
calculations gave S = 2.38 0, 6 = 2.7 0,(1,0, S < 2.01 
7 1 .while finite- T simulations obtain 5 between 2 and 4.8 
lEIi from the filling of the gap as g(fi) cx T s @, ©, OH . 

It was also suggested long ago jTl| that disordered in- 
sulators enter a glass state at low temperature. Ample 
experimental and numerical evidence of glassy nonequi- 
librium effects in these systems has been obtained since 
[l2l ]. However, it remains unclear whether these effects 
are purely dynamical or reflect an underlying transition 
to an equilibrium glass phase (CP), and whether there is 
a link between glassiness and the Coulomb gap. Some 
evidence for a sharp equilibrium transition to a GP was 
found in simulations of localized charges with random 
positions liL fll . fl5l | but not in the presence of on-site 
disorder [7J, |l5(. In the latter case, the transition would 
not break any symmetry of the Hamiltonian, similarly 
to the l ong -debated Almeida-Thouless transition in spin 
glasses [l6( . These issues have been brought a gain to the 
fore by recent mean-field studies 13, 17, Il8l Il9l| which 
predict a "replica symmetry broken" equilibrium GP be- 
low a critical temperature T g in the presence of on-site 
disorder. In this GP correlations remain critical, which 
leads to S = d — 1, and both T g and the gap width A 
scale as W~? for d — 3 and large disorder strength W 
0. 




FIG. 1: (Color online) Phase diagram of the Coulomb glass 
model. The thin lines are the simulation paths [BC: T = 
(3/10)W + 7/100, Wc = 0.977; DE: T = (9/80) W + 3/200, 
We = 0.506; T A = 0.275.] The fluid-COP boundary inter- 
polates the transition temperatures estimated along AB, BC, 
DE, shown in red. Our results indicate that no glass phase 
exists above the blue dashed line. See also Fig. 2 of Ref . [l8| . 



In this Letter, we investigate these predictions via ex- 
tensive Monte Carlo (MC) simulations of the Coulomb 
glass lattice model with on-site disorder [2(|. In addi- 
tion, we study in detail the transition from the fluid to 
the charge-ordered phase (COP). For W = 0, there is 
good numerical evidence for an Ising-like transition [2l| . 
For W 7^ 0, mean- field theory predicts a stable COP for 
d = 3 [18j,|22|. Beyond mean field there is some numerical 
evidence that the COP survives small positional disorder 
d, 0] and on-site disorder [23| , but neither the phase di- 
agram nor the critical properties have been investigated. 
Our results are as follows: (i) A COP exists below the 
(approximate) phase boundary in Fig. 1. (ii) The fluid- 
COP transition is consistent with the Random Field Ising 
model (RFIM) universality class, which shows that the 
interaction is effectively screened near the transition, (iii) 
No GP is found for T well below the mean- field T g [18j |. 
in agreement with the results of Ref. 0] but at lower T 
and in a wider range for W. (iv) Due to the long-range 
interaction, the glass correlation length increases rapidly 
and possibly diverges as T — > 0. (v) For large W, the 
DOS scales as g(e, T) = T 5 f(\e\/T) near the gap, with a 
saturated exponent S ~ 2.00. (vi) As W decreases, scal- 
ing breaks down above the COP, and an effective power 
law <?l(£, T = 0) cx \e\ s holds with 5 > 3, in contrast with 
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Ref. A more extended account will appear later [2! 
Model and simulation - We study the Hamiltonian 
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-^fa-JQ H-W^n,^ (1) 



where ni G {0, 1} are the occupation numbers for the 
N = L d sites of a hypercubic lattice (d = 3) with 
Sili n i = KN, and is the distance from i to j. The 
filling factor is K = 1/2 (which gives /u = 0). The ran- 
dom on-site energies ipi are independen and Gaussian- 
distributed with zero mean and variance unity. Energies 
and temperatures will be in units of e 2 / (n£) and lengths 
in units of the lattice spacing £. 

We carry out canonical MC sampling along the paths 
ABC and ADE in Fig. 1 and at constant W — 
0.2, 0.5, 1, 2, 4. We consider an infinite sphere of periodic 
images of a central L 3 cell and sum over all interactions 
with the Ewald method with a dipole surface term [25j |. 
To reach low temperatures, we use the exchange MC algo- 
rithm |26] . For each realization (sample) ip = {(pi}jLi, we 
simulate identical replicas with different (T, W) along the 
simulation path. Every N/2 Metropolis steps for single- 
electron hops, replicas at adjacent (T, W), (T", W) swap 
their configurations with probability min(l,p), where p = 
exp[(/3 - p')(H - W) + {W - W)(pR! - p"R)], f3 = 1/T, 
and 1Z = J2i n iVi' which preserves detailed balance. The 
simulation time t s is chosen so that averages over the 
intervals [t s /3,t s ] and [t s /9, t s /3] agree within the sta- 
tistical errors, and that the identity 2TN~ 1 {(lZ)] av = 
W(2N-^l 1 [(n^nf ) )] av -l 
order, is satisfied. Here, (•} and 
sample averages and a, b are two independently simulated 
replicas with the same (<p,T, W) [27l | . 

Charge ordering - Fig. 2 (top inset) shows the COP 
order parameter M s — [{\m s \)] av along the paths AB, 
BC, and DE, where m s — N^ 1 J^ i=1 o~i and o~i — 
Si(—l) Xi+Vi+Zi (we introduce the Ising variables Sj = 
2ni — 1). The sharp increase demonstrates a transition 
to a COP. To determine the transition temperature T c , 
we measure the finite-size correlation length (CL) (28| 
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where Xi( k ) = N [ J2i j[(°'i a 'i)]a^e ik ' r « and k min = 
(2n/L, 0, 0). Along BC, the data for £ L (T) /L for different 
L cross (Fig. 2, main panel), which signals [2§| a transi- 
tion at T C BC = 0.0950(15). We observe similar crossings 
along AB and DE (not shown) at T^ B = 0.1280(15), 
in excellent agreement with Refs.(6l. \2~]\, and T C DE = 
0.031(2). The curve (T c (T^)/T c AB ) l fi0 = l-(W/0.15) l eo 
interpolates these three points and gives the approximate 
fluid-COP phase boundary in Fig. 1. 

Critical behavior - Since at W = the fluid-COP tran- 
sition has a positive specific-heat exponent [2l( , disorder 




FIG. 2: (Color online) Charge-order CL along path BC in 
Fig. 1. Top inset: order parameter M s along paths AB, BC, 
and DE. Bottom inset: specific heat along path BC. 

TABLE I: Critical exponents for the fluid-COP transition 
along BC in Fig. 1, compared with the RFIM values [32| . 
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Coulomb glass 1.69(17) 


2.89(9) 


0.06(4) 


1.11(12) 


RFIM 1.44(12) 


2.93(11) 


0.011(4) 


1.37(9) 



is relevant and the W ^ transition will be governed by 
a random fixed point which, by analogy with the RFIM 
[2ll | . we expect to be at T = [3(|. Assuming that the 
W ^ transition is second order (indeed the distribu- 
tion of m s is unimodal at all T for a predominant, and 
increasing with L, fraction of the samples [24]]) we ob- 
tain the critical exponents in Table 1.3/u and *y/v were 
estimated with the quotient method [31| for the observ- 
ables M s and Xl — N[(m s 2 }) av respectively [the quo- 
tient estimates from {L,L r ) = (6, 8), (6, 10) and (8,10) 
agree within the errors], while 7/f was obtained by fit- 
ting aU 1 l v to the height of the peak of the susceptibility 



N[(m, s 2 ) — (\m s \) 2 ]av (data not shown). The peak height 
for the specific heat c L = 1/{NT 2 )[(H 2 ) - {H) 2 ] av in- 
creases slowly with L (Fig. 2, bottom inset), which sug- 
gests either a < or a logarithmic divergence (a = 0). 
We could not directly estimate v in a reliable way, but 
we obtain v = 1.11(12) from the modified hyperscaling 
relation (30j (d — ff)v = 2 — a, assuming a — and us- 
ing 9 = j/v - 7/1/ = 1.20(20). As shown in Table I, 
the exponents ag ree fairly well with the known values 
for the RFIM [32| , which suggests that the interaction is 
effectively short-range near the phase boundary. 

Glass phase - Several works have searched for a GP by 
measuring the parameter [((rii) — l/2) 2 ] av ll| or higher 
cumulants of the overlap between two replicas 15(. We 
measure instead the glass CL £® obtained from Eq. (2) 
by replacing [(<Xj<7j)] at , with the "spin-glass" correlation 
function G(r i:j ) = [{(SiSj) - (Si)(Sj)) 2 ] av . In the fluid 
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phase we have G(r) ~ exp(— r/£ G ) for £ G < r <C L, 
where £ G is the bulk CL, thus ~ £ G for L > £ G 
and £1 ~ i for L <C £ G . In a "many-state" GP [la ]. 
G(r) tends to a constant for large r, thus we have £° ~ 
Hence the existence of a GP will be signaled 
by the crossing of ££(T)/£ for different L near T = T g 
[29(. As shown in Fig. 3(a) for W — 0.5, 1, we observe no 




FIG. 3: (Color online) (a) Glass correlation length £ G at 
W = 1 and W = 0.5 (shifted upwards by a factor 7). The 
absence of crossing is evidence against the existence of an 
equilibrium glass transition, (b) Power-law fit of £ G o = cT~" 
for W = 1. (c) Scaling plot = Lf(TL 1/v ') for W = 1 and 
W = 0.5 (shifted to the right by a factor 5). 

crossing down to the lowest equilibrated temperature and 
well below the mean-held glass transition (T g rj 0.037 for 
W = 0.5 [HI). Similar results were found in Ref.0] for 
T > 0.03 and W < 0.4. We also exclude that a GP 
occurs at T g > T c along BC and DE by comparing the 
crossing temperatures for £,1/L (not shown) and £l/L: 
For all pairs (L,L') they differ by less than 1%. Together 
with the results at constant W , this indicates that no GP 
exists above the dashed line in Fig. 1. 

Fig. 3(b) shows that is nearly independent of L at 
large T and W = 1 apart from small finite-size effects, 
thus Q =w is a good estimate of £ G for T > 0.017. At 
lower T we observe ~ L, which shows that £ G be- 
comes larger than ~ 10, and possibly diverges as T — > 0. 
Indeed, £° =10 (T) can be fitted for T £ [0.017,0.113] by 
a power law T~ u ' with v 1 = 0.82 [Fig. 3(b)], which also 
gives a satisfactory finite-size scaling £° = Lj(TL 1 l v ) 
[Fig. 3(c)] . A divergence would be a nontrivial prediction, 
since for W ^ the ground state of a periodic sample is 
unique. Because L and are rather small, however, we 
cannot rule out neither that £ G stays finite at T — 0, nor 
an exponent v 1 = 1. An interesting question is whether 
there is a link with the T~ l divergence of the screening 
length found in mean-field theory [l3T ] and from simple 



arguments H3, El- We tested that the CL obtained by 
replacing [{(Ti(Tj)] av with [((TiO-j) - (?i){vj)]av in Eq. (2), 
remains smaller than unity at all T, which suggests that 
the correlated regions are disordered. Finally, we simu- 
lated the short-range RFIM choosing W so that the value 
of N~ 1 {{J2 i TH<Pi)]av is close to the Coulomb glass value 
at W — 1, and found < 1 as T — > 0, which suggests 
that the large CL is due to the long-range interaction. 

Coulomb gap - The single-particle DOS is defined as 
g L {e,T) =iV" 1 [(EiLi<K e - e i)}W where e t = E,-^ - 



K)/t 



Wifi is the cost of adding an electron at site i 



of the central cell while leaving the periodic images un- 
changed. We compute the infinite sum with the Ewald 
method. Because of the dipole term 25j, the DOS has 
no hard gap 24j, |34[, unlike for a finite, nonperiodic sys- 
tem. The finite-size effects due to the energy scale L" 1 
turn out to be significant for |e| < aL^ 1 with a ~ 0.3, 
while those due to the sample fluctuations of /i (of order 
W/L d / 2 ) were drastically reduced by shifting the DOS 
before averaging over the samples [351 ] . In the gap region 
(|e|, T) <C A, which is our only focus here, one expects the 
scaling g L {e,T) = T 5 f{\e\/T) for ai" 1 < (|e|,T), with 
f(x) ~ constant as x — > and f(x) ~ ex 5 as x — > oo. 
Fig. 4(a,b,c) show scaling plots with <5 = 2 for L = 10 
and W — 4,2, 0.5. For W — 4 scaling is excellent even for 
this moderate size, with small deviations for |e| < 0.03 
due to finite-size effects. For A/T^|e|/T>6 the data 
are well fitted by gio(e, T) = c|e| 2 with c ~ 1.1, which is 
close to the self-consistent prediction c = 3/tt [20] (while 
Ref.[lJ] finds c = 0.2083). As shown in Fig. 4(b) (inset), 
the finite-size scaling ansatz t/i(e,T) = L~ s h(eL), which 
should hold for T < |e| < aL^ 1 < A [with h(x) ~ c\x\ s 
for large x], is also well satisfied with c = 3/tt, 5 = 2. Our 
final estimate is 5 — 2.01 ± 0.05, which provides strong 
support for a saturated ES bound. 

For decreasing W, we observe increasingly stronger de- 
viations from the 8 — 2 scaling. A fit gio(e,T) = c\e\ s at 
low T gives an effective exponent S ~ 2.3 for W = 2 
and S > 2.8 for W = 0.5 [Fig. 4(c)]. We interpret 
this as a crossover due to the vicinity of the fluid-COP 
boundary, below which the DOS has a hard gap at 
T = 0. The crossover is apparent in Fig. 4(d): Since 
9L ( e = 0,T)/T 2 cx T 8 - 2 for aL" 1 < T < A, the plateau 
for W = 4 supports 5 = 2, while for decreasing W the 
exponent increases to S > 3. The L dependence is con- 
sistent with the scaling g L (e = 0,T) = T s h(TL) (not 
shown) with 5 extracted from Fig. 4(d) for each value of 
W. Our results differ markedly from Ref.Q, which re- 
ports S = 1.83(3) for the same model at W = 0.4, T = 
[however, <?l(£ — 0) is much larger than our data, and 
increases with L], A similar crossover in the DOS was 
reported for d — 2, where the COP occurs at W = [3t3 ]. 

In conclusion, we presented evidence that no equilib- 
rium glass phase exists in the Coulomb glass, but a sat- 
urated ES bound holds. The long-range part of the in- 
teraction appears to be irrelevant as to the equilibrium 
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FIG. 4: (Color online) (a,b,c) Scaling plots of the DOS 
for W = 4, 2, 0.5 and L = 10. From bottom to top: T = 
0.047, 0.035, 0.025, 0.017, 0.0105 plus T = 0.0049 for W = 4, 2 
and T = 0.0026 for W = 4. The error is the standard devia- 
tion of the sample fluctuations. The solid lines represent the 
ES law g(e, T) = 3e 2 /%. The dashed line with slope 2.8 in 
(c) highlights the departure from the ES law at low W (the 
W — 0.5, T = 0.0105 data are not fully equilibrated, but the 
slope increases with simulation time). Inset of (b) Finite- 
size scaling for W = 4 and T = 0.0077 (other values of W 
and T give similar plots), (d) Temperature dependence of 
the DOS at |e| < 0.0075 for L = 10. 

thermodynamics, except for a possible diverging correla- 
tion length at T = 0, which calls for further investigation. 
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